function TS = volTS(R, HL, annual_factor)
% volTS: Estimate time varying volatilities
% TS = volTS(R, halflife, annual_factor)  Estimate time varying volatilities
% using simple expanding window if no halflife provided.  EWMA otherwise.
% Example:
%    TS = volTS(R, halflife, annual_factor);
%
% TODO: add de-meaning functionality
%

N = length(R.dates);


% Build initial vol using halflife
if exist('HL','var')
    gamma = 0.5^(1/HL);
    decaypwr = (HL-1:-1:0);
    wgts = (gamma).^decaypwr;
    wgts = wgts/sum(wgts);
    init_vol = wgts*(R.data(1:HL,:).^2);
else
    error('no functionality for missing halflife yet.');
end

% Build covariance matrix time series
TS = buildTS([], R.header, R.dates);  % volatility TS structure

last_vol_t = init_vol;
for t = 1 : N

   rt = R.data(t, :); % Single period returns

   % Estimate volatilities through time
   vol_t = gamma * last_vol_t + (1 - gamma) * (rt .* rt);

   % Update t-1 values
   last_vol_t = vol_t;
   
   % Save current time period outputs
   TS.data(t, :) = sqrt(vol_t);
   
end

% Annualize
if exist('annual_factor','var')
    TS = multTS(TS, sqrt(annual_factor));
end
% Lag cov matrix for use in backtesting
TS.data = cat(1, TS.data(1, :), TS.data(1:end-1,:));

return
